home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / mkconfig.c < prev    next >
C/C++ Source or Header  |  1988-11-24  |  11KB  |  360 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: mkconfig.c,v 1.4 85/08/26 10:41:39 timo Exp $
  5. */
  6.  
  7. /* Generate constants for configuration file            */
  8.  
  9. /* If your C system is not unix but does have signal/setjmp,    */
  10. /*    add a #define unix                    */
  11. /* You may also need to add some calls to signal().        */
  12.  
  13. #ifdef unix
  14.  
  15. #define SIGNAL
  16.  
  17. #include <signal.h>
  18. #include <setjmp.h>
  19.  
  20.     jmp_buf lab;
  21.     overflow(sig) int sig; { /* what to do on overflow/underflow */
  22.         signal(sig, overflow);
  23.         longjmp(lab, 1);
  24.     }
  25.  
  26. #else
  27.     /* Dummy routines instead */
  28.     int lab=1;
  29.     int setjmp(lab) int lab; { return(0); }
  30.  
  31. #endif
  32.  
  33. #define fabs(x) (((x)<0.0)?(-x):(x))
  34. #define min(x,y) (((x)<(y))?(x):(y))
  35.  
  36. /* These routines are intended to defeat any attempt at optimisation */
  37. Dstore(a, b) double a, *b; { *b=a; }
  38. double Dsum(a, b) double a, b; { double r; Dstore(a+b, &r); return (r); }
  39. double Ddiff(a, b) double a, b; { double r; Dstore(a-b, &r); return (r); }
  40. double Dmul(a, b) double a, b; { double r; Dstore(a*b, &r); return (r); }
  41. double Ddiv(a, b) double a, b; { double r; Dstore(a/b, &r); return (r); }
  42.  
  43. double power(x, n) int x, n; {
  44.     double r=1.0;
  45.     for (;n>0; n--) r*=x;
  46.     return r;
  47. }
  48.  
  49. int log(base, x) int base; double x; {
  50.     int r=0;
  51.     while (x>=base) { r++; x/=base; }
  52.     return r;
  53. }
  54.  
  55. main(argc, argv) int argc; char *argv[]; {
  56.     char c;
  57.     short newshort, maxshort, maxershort;
  58.     int newint, maxint, shortbits, bits, mantbits,
  59.         *p, shortpower, intpower, longpower;
  60.     long newlong, maxlong, count;
  61.  
  62.     int i, ibase, iexp, irnd, imant, iz, k, machep, maxexp, minexp,
  63.         mx, negeps, ngrd, normalised;
  64.     double a, b, base, basein, basem1, eps, epsneg, xmax, newxmax,
  65.            xmin, xminner, y, y1, z, z1, z2;
  66.  
  67.     double BIG, Maxreal;
  68.     int BASE, MAXNUMDIG, ipower, tenlogBASE, Maxexpo, Minexpo, Dblbits;
  69.  
  70. #ifdef SIGNAL
  71.     signal(SIGFPE, overflow);
  72.     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
  73. #endif
  74.  
  75. /****** Calculate max short *********************************************/
  76. /*      Calculate 2**n-1 until overflow - then use the previous value    */
  77.  
  78.     newshort=1; maxshort=0;
  79.  
  80.     if (setjmp(lab)==0)
  81.         for(shortpower=0; newshort>maxshort; shortpower++) {
  82.             maxshort=newshort;
  83.             newshort=newshort*2+1;
  84.         }
  85.  
  86.     /* Now for those daft Cybers: */
  87.  
  88.     maxershort=0; newshort=maxshort;
  89.  
  90.     if (setjmp(lab)==0)
  91.         for(shortbits=shortpower; newshort>maxershort; shortbits++) {
  92.             maxershort=newshort;
  93.             newshort=newshort+newshort+1;
  94.         }
  95.  
  96.     bits= (shortbits+1)/sizeof(short);
  97.     c= (char)(-1);
  98.     printf("/\* char=%d bits, %ssigned *\/\n", sizeof(c)*bits,
  99.             ((int)c)<0?"":"un");
  100.     printf("/\* maxshort=%d (=2**%d-1) *\/\n", maxshort, shortpower);
  101.  
  102.     if (maxershort>maxshort) {
  103.         printf("/\* There is a larger maxshort, %d (=2**%d-1), %s *\/\n",
  104.             maxershort, shortbits, 
  105.             "but only for addition, not multiplication");
  106.     }
  107.  
  108. /****** Calculate max int by the same method ***************************/
  109.  
  110.     newint=1; maxint=0;
  111.  
  112.     if (setjmp(lab)==0)
  113.         for(intpower=0; newint>maxint; intpower++) {
  114.             maxint=newint;
  115.             newint=newint*2+1;
  116.         }
  117.  
  118.     printf("/\* maxint=%d (=2**%d-1) *\/\n", maxint, intpower);
  119.  
  120. /****** Calculate max long by the same method ***************************/
  121.  
  122.     newlong=1; maxlong=0;
  123.  
  124.     if (setjmp(lab)==0)
  125.         for(longpower=0; newlong>maxlong; longpower++) {
  126.             maxlong=newlong;
  127.             newlong=newlong*2+1;
  128.         }
  129.  
  130.     if (setjmp(lab)!=0) { printf("\nUnexpected under/overflow\n"); exit(1); }
  131.  
  132.     printf("/\* maxlong=%ld (=2**%d-1) *\/\n", maxlong, longpower);
  133.  
  134. /****** Pointers ********************************************************/
  135.     printf("/\* pointers=%d bits%s *\/\n", sizeof(p)*bits,
  136.         sizeof(p)>sizeof(int)?" BEWARE! larger than int!":"");
  137.  
  138. /****** Base and size of mantissa ***************************************/
  139.     a=1.0;
  140.     do { a=Dsum(a, a); } while (Ddiff(Ddiff(Dsum(a, 1.0), a), 1.0) == 0.0);
  141.     b=1.0;
  142.     do { b=Dsum(b, b); } while ((base=Ddiff(Dsum(a, b), a)) == 0.0);
  143.     ibase=base;
  144.     printf("/\* base=%d *\/\n", ibase);
  145.  
  146.     imant=0; b=1.0;
  147.     do { imant++; b=Dmul(b, base); }
  148.     while (Ddiff(Ddiff(Dsum(b,1.0),b),1.0) == 0.0);
  149.     printf("/\* Significant base digits=%d *\/\n", imant);
  150.  
  151. /****** Various flavours of epsilon *************************************/
  152.     basem1=Ddiff(base,1.0);
  153.     if (Ddiff(Dsum(a, basem1), a) != 0.0) irnd=1; 
  154.     else irnd=0;
  155.  
  156.     negeps=imant+imant;
  157.     basein=1.0/base;
  158.     a=1.0;
  159.     for(i=1; i<=negeps; i++) a*=basein;
  160.  
  161.     b=a;
  162.     while (Ddiff(Ddiff(1.0, a), 1.0) == 0.0) {
  163.         a*=base;
  164.         negeps--;
  165.     }
  166.     negeps= -negeps;
  167.     printf("/\* Smallest x such that 1.0-base**x != 1.0=%d *\/\n", negeps);
  168.  
  169.     epsneg=a;
  170.     if ((ibase!=2) && (irnd==1)) {
  171.     /*    a=(a*(1.0+a))/(1.0+1.0); => */
  172.         a=Ddiv(Dmul(a, Dsum(1.0, a)), Dsum(1.0, 1.0));
  173.     /*    if ((1.0-a)-1.0 != 0.0) epsneg=a; => */
  174.         if (Ddiff(Ddiff(1.0, a), 1.0) != 0.0) epsneg=a;
  175.     }
  176.     printf("/\* Small x such that 1.0-x != 1.0=%g *\/\n", epsneg);
  177.     /* it may not be the smallest */
  178.  
  179.     machep= -imant-imant;
  180.     a=b;
  181.     while (Ddiff(Dsum(1.0, a), 1.0) == 0.0) { a*=base; machep++; }
  182.     printf("/\* Smallest x such that 1.0+base**x != 1.0=%d *\/\n", machep);
  183.  
  184.     eps=a;
  185.     if ((ibase!=2) && (irnd==1)) {
  186.     /*    a=(a*(1.0+a))/(1.0+1.0); => */
  187.         a=Ddiv(Dmul(a, Dsum(1.0, a)), Dsum(1.0, 1.0));
  188.     /*    if ((1.0+a)-1.0 != 0.0) eps=a; => */
  189.         if (Ddiff(Dsum(1.0, a), 1.0) != 0.0) eps=a;
  190.     }
  191.     printf("/\* Smallest x such that 1.0+x != 1.0=%g *\/\n", eps);
  192.  
  193. /****** Round or chop ***************************************************/
  194.     ngrd=0;
  195.     if (irnd == 1) { printf("/\* Arithmetic rounds *\/\n"); }
  196.     else { 
  197.         printf("/\* Arithmetic chops");
  198.         if (Ddiff(Dmul(Dsum(1.0,eps),1.0),1.0) != 0.0) {
  199.             ngrd=1;
  200.             printf(" but uses guard digits");
  201.         }
  202.         printf(" *\/\n");
  203.     }
  204.  
  205. /****** Size of and minimum normalised exponent ****************************/
  206.     y=0; i=0; k=1; z=basein; z1=(1.0+eps)/base;
  207.  
  208.     /* Coarse search for the largest power of two */
  209.     if (setjmp(lab)==0) /* in case of underflow trap */
  210.         do {
  211.             y=z; y1=z1;
  212.             z=Dmul(y,y); z1=Dmul(z1, y);
  213.             a=Dmul(z,1.0);
  214.             z2=Ddiv(z1,y);
  215.             if (z2 != y1) break;
  216.             if ((Dsum(a,a) == 0.0) || (fabs(z) >= y)) break;
  217.             i++;
  218.             k+=k;
  219.         } while(1);
  220.  
  221.     if (ibase != 10) {
  222.         iexp=i+1; /* for the sign */
  223.         mx=k+k;
  224.     } else {
  225.         iexp=2;
  226.         iz=ibase;
  227.         while (k >= iz) { iz*=ibase; iexp++; }
  228.         mx=iz+iz-1;
  229.     }
  230.  
  231.     /* Fine tune starting with y and y1 */
  232.     if (setjmp(lab)==0) /* in case of underflow trap */
  233.         do {
  234.             xmin=y; z1=y1;
  235.             y=Ddiv(y,base); y1=Ddiv(y1,base);
  236.             a=Dmul(y,1.0);
  237.             z2=Dmul(y1,base);
  238.             if (z2 != z1) break;
  239.             if ((Dsum(a,a) == 0.0) || (fabs(y) >= xmin)) break;
  240.             k++;
  241.         } while (1);
  242.  
  243.     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
  244.  
  245.     minexp= (-k)+1;
  246.  
  247.     if ((mx <= k+k-3) && (ibase != 10)) { mx+=mx; iexp+=1; }
  248.     printf("/\* Number of bits used for exponent=%d *\/\n", iexp);
  249.     printf("/\* Minimum normalised exponent=%d *\/\n", minexp);
  250.     printf("/\* Minimum normalised positive number=%g *\/\n", xmin);
  251.  
  252. /****** Minimum exponent ***************************************************/
  253.     if (setjmp(lab)==0) /* in case of underflow trap */
  254.         do {
  255.             xminner=y;
  256.             y=Ddiv(y,base);
  257.             a=Dmul(y,1.0);
  258.             if ((Dsum(a,a) == 0.0) || (fabs(y) >= xminner)) break;
  259.         } while (1);
  260.  
  261.     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
  262.  
  263.     normalised=1;
  264.     if (xminner != 0.0 && xminner != xmin) {
  265.         normalised=0;
  266.         printf("/\* The smallest numbers are not kept normalised *\/\n");
  267.         printf("/\* Smallest unnormalised positive number=%g *\/\n",
  268.             xminner);
  269.     }
  270.  
  271. /****** Maximum exponent ***************************************************/
  272.     maxexp=2; xmax=1.0; newxmax=base+1.0;
  273.     if (setjmp(lab) == 0) {
  274.         while (xmax<newxmax) {
  275.             xmax=newxmax;
  276.             newxmax=Dmul(newxmax, base);
  277.             if (Ddiv(newxmax, base) != xmax) break; /* ieee infinity */
  278.             maxexp++;
  279.         }
  280.     }
  281.     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
  282.  
  283.     printf("/\* Maximum exponent=%d *\/\n", maxexp);
  284.  
  285. /****** Largest and smallest numbers ************************************/
  286.     xmax=Ddiff(1.0, epsneg);
  287.     if (Dmul(xmax,1.0) != xmax) xmax=Ddiff(1.0, Dmul(base,epsneg));
  288.     for (i=1; i<=maxexp; i++) xmax=Dmul(xmax, base);
  289.     printf("/\* Maximum number=%g *\/\n", xmax);
  290.  
  291. /****** Hidden bit + sanity check ***************************************/
  292.     if (ibase != 10) {
  293.         mantbits=log(2, (double)ibase)*imant;
  294.         if (mantbits+iexp+1 == sizeof(double)*bits+1) {
  295.             printf("/\* Double arithmetic uses a hidden bit *\/\n");
  296.         } else if (mantbits+iexp+1 == sizeof(double)*bits) {
  297.             printf("/\* Double arithmetic doesn't use a hidden bit *\/\n");
  298.         } else {
  299.             printf("/\* Something fishy here! %s %s *\/\n",
  300.                 "Exponent size + mantissa size doesn't match",
  301.                 "with the size of a double.");
  302.         }
  303.     }
  304.  
  305. /****** The point of it all: ********************************************/
  306.     printf("\n/\* Numeric package constants *\/\n");
  307.  
  308.     BIG= power(ibase, imant)-1.0;
  309.     MAXNUMDIG= log(10, BIG);
  310.     ipower= log(10, (double)(maxint/2));
  311.     Maxreal= xmax;
  312.     Maxexpo= log(2, (double)ibase)*maxexp;
  313.     Minexpo= log(2, (double)ibase)*minexp;
  314.     Dblbits= log(2, (double)ibase)*imant;
  315.     tenlogBASE= min(MAXNUMDIG/2, ipower);
  316.     BASE=1; for(i=1; i<=tenlogBASE; i++) BASE*=10;
  317.  
  318.     printf("#define BIG %1.1f /\* Maximum integral double *\/\n", BIG);
  319.     printf("#define LONG ");
  320.         for(i=1; i<=MAXNUMDIG; i++) printf("9");
  321.         printf(".5 /\* Maximum power of ten less than BIG *\/\n");
  322.     printf("#define MAXNUMDIG %d /\* The number of 9's in LONG *\/\n",
  323.         MAXNUMDIG);
  324.     printf("#define MINNUMDIG 6 /\* Don't change: this is here for consistency *\/\n");
  325.  
  326.     printf("#define BASE %d /\* %s *\/\n",
  327.         BASE, "BASE**2<=BIG && BASE*2<=Maxint && -2*BASE>=Minint");
  328.     if (((double)BASE)*BASE > BIG || ((double)BASE)+BASE > maxint) {
  329.         printf("*** BASE value wrong\n");
  330.         exit(1);
  331.     }
  332.     printf("#define tenlogBASE %d /\*  = log10(BASE) *\/\n", tenlogBASE);
  333.     printf("#define Maxreal %e /\* Maximum double *\/\n", Maxreal);
  334.     printf("#define Maxexpo %d /\* Maximum value such that 2**Maxexpo<=Maxreal *\/\n",
  335.         Maxexpo);
  336.     printf("#define Minexpo (%d) /\* Minimum value such that -2**Minexpo>=Minreal *\/\n",
  337.         Minexpo);
  338.     printf("#define Dblbits %d /\* Maximum value such that 2**Dblbits<=BIG *\/\n",
  339.         Dblbits);
  340.  
  341.     printf("#define Maxintlet %d /\* Maximum short *\/\n", maxshort);
  342.     printf("#define Maxint %d /\* Maximum int *\/\n", maxint);
  343.  
  344.     printf("#define Maxsmallint %d /\* BASE-1 *\/\n", BASE-1);
  345.     printf("#define Minsmallint (%d) /\* -BASE *\/\n", -BASE);
  346.  
  347. /* An extra goody: the approximate amount of data-space */
  348. /* Put here because it is likely to be slower then the rest */
  349.  
  350.     /*Allocate blocks of 1000 until no more available*/
  351.     /*Don't be tempted to change this to 1024: */
  352.     /*we don't know how much header information there is*/
  353.  
  354.     for(count=0; (p=(int *)malloc(1000))!=0; count++) { }
  355.  
  356.     printf("\n/\* Memory~= %d000 *\/\n", count);
  357.  
  358.     exit(0);
  359. }
  360.